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ABSTRACT 

. Context. Numerical simulations of planets embedded in protoplanetary gaseous discs are a precious tool for studying the planetary 

migration ; however, some approximations have to be made. Most often, the selfgravity of the gas is neglected. In that case, it is not 
clear in the literature how the material inside the Roche lobe of the planet should be taken into account. 

Aims. Here, we want to address this issue by studying the influence of various methods so far used by different authors on the 
Q . migration rate. 

'_ , ■ Methods. We performed high-resolution numerical simulations of giant planets embedded in discs. We compared the migration rates 

' with and without gas selfgravity, testing various ways of taking the circum-planetary disc (CPD) into account. 

Results. Different methods lead to significantly different migration rates. Adding the mass of the CPD to the perturbing mass of 
the planet accelerates the migration. Excluding a part of the Hill sphere is a very touchy parameter that may lead to an artificial 
suppression of the type III, runaway migration. In fact, the CPD is smaller than the Hill sphere. We recommend excluding no more 
than a 0.6 Hill radius and using a smooth filter. Alternatively, the CPD can be given the acceleration felt by the planet from the rest of 
the protoplanetary disc. 

Conclusions. The gas inside the Roche lobe of the planet should be very carefully taken into account in numerical simulations without 
any selfgravity of the gas. The entire Hill sphere should not be excluded. The method used should be explicitly given. However, no 
method is equivalent to computing the full selfgravity of the gas. 

O ■ 

' Key words. Methods: numerical — (Stars): planetary systems: formation — Accretion, accretion discs 

o 

o\ . 

f^} 1. Introduction an overall consistency between the codes; however, a few nu- 

., , ,. , • , , , merical and physical issues are still open. In particular, a circum- 
Planetary migration has been widely studied in the past decade, lanet disc (CPD) generally appears around giant planets. 
^ ; after the discovery of the first exoplanets. Indeed, the first known without mesh refinement; this stmcture is generally poorly re- 
, exoplanets were mostly hot Jupiters that could not have formed solyed> but m t Afj k ^ ^ tQ ^ ^ h ^ 
Oj ■ where they presently orbit, according to the standard theory of tentiaU eX ert a strong force on it. So far, a very limited num- 
planet formation. Consequently, they are supposed to have mi- ber of smdies have considered the probably important role of the 
grated in the protoplanetary gaseous disc before its dissipation. CpD fa ^ mi tion rate As this material is consi dered to be- 
More recently, pairs of exoplanets in mean motion resonance , tQ ^ ^ and ^ order tQ ayoid numerical due to 
have been detected. This particular configuration requires some ^ ^ ^ SQme authors exdude ^ m] here of ^ 
dissipative process for the planets to approach each other and lanet from ^ cakulation of the force of the on the lanet 
lock in resonance. This is a strong evidence for planetary migra- Mogt authors identify ±e QpD wkh ^ ffil] s phere , which in 

tlon ' , . . fact does not hold as we shall see below (see also D' Ange lo et al.l 

The theory of planetary migratio n was in fact pro- j^k Some authors add the mass of the - m the ffin here tQ 

£9ied before the observa t ional e vidence (| G o ldi-eich & Tremaind ^ lanet mass whgn co tin the plane t gravity field. Others 
1979, Lin & Emato UmSyto^aa & Sicardy| | l987t add the acce i er ation caused by the entire gas disc on the planet 
Lin & Papaloizou| | l986bt Ward| | 1986|) . Numerical simulations tQ ^ ^ ^ m] here ^ f ^ skuation ^ confused ^ 



allow study of more complex cases (several planets for instance) and there ^ nQ consensus on this issue, 
and refinement of our understanding of this phenomenon. 

Therefore, numerica l studies of planetar y migration are numer- , . . . 

ous in the literature. iDe Val-Borro etaTI (120061) have compared Recently, [Crida etalj d2008|) have shown how planetary mi- 

the results of various codes on two typical test problems and find g ratlon can ex P lain the orbltal Parameters of exoplanets in mean 

motion resonance, and remarked that the outcome of their sim- 

Send offprint requests to: A. Crida, ulations had a dependence on the way they deal with the gas in 

crida@tat.physik.uni-tuebingen.de the Hill sphere of the planets. Pepliri ski et al.l d2008l) have also 
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enlightened that the role of the CPD is not a minor one, and that 
a careful study of this problem is needed to answer this question. 

In this paper, we study the influence of the CPD on the plane- 
tary migration in numerical simulations. The theory is discussed, 
to get a better understanding of the problem. This leads to pos- 
sible solutions to treat the CPD. We perform 2D simulations 
with modest to high resolution, and we test and compare var- 
ious recipes. Our goal is to provide a survey of the parameter 
"CPD" in numerical simulations. In particular, the fraction of 
the Hill sphere that is excluded is a very touchy parameter, and 
the way the exclusion is performed has a great influence on the 
outcome of the simulation. Many authors are not aware of this. 
Our survey shows that some popular methods are inappropriate. 
We present better ways of dealing with the material present in 
the Hill sphere, in particular the CPD. 

This requires an analysis of the Hill sphere structure, but in 
the framework of global disc simulations, which are not dedi- 
cated to resolve the Hill sphere with a high level of precision. 
We do not aim at understanding the detailed physics and the fine 
structure of the CPD itself, but more how the gas in the Hill 
sphere of the planet should be taken into account. Convergence 
properties and the impact of varied prescriptions on the equation 
of state (EOS) or the softening length are fundamental issues, but 
we do not address them in this work. Instead, our aim is to find 
an adequate prescription for the torque calculation in non self- 
gravitating simulations, based on general considerations such as 
the conservation of angular momentum and the main features of 
the flow topology in the vicinity of the planet. The conservation 
of angular momentum is enforced in our scheme independently 
of the resolution, while the flow topology depends very weakly 
on the resolution (even if the mass of the CPD does, in an EOS- 
dependent way). 

This paper is organised as follows. In Sect. [2] the problem is 
explained and analysed. Some definitions are provided, in order 
to clarify the strategy and the questions that arise from the pres- 
ence of a CPD. We show that the problem can be divided into the 
questions of the perturbing mass, of the gravitational mass, and 
of the inertial mass of the planet. In Sect. [3] the numerical setup 
is presented. Then, in the following sections, numerical experi- 
ments are performed. In Sect.|H a Jupiter mass planet, in type II 
migration in considered. In Sect. [5] the influence of the CPD on 
a Saturn mass planet, in type III migration is studied. In Sect. [6] 
the structure of the Hill sphere is analysed in detail, and the size 
of the CPD is estimated. Last, our results are summarised, and 
we conclude in Sect. [7] 

2. Problem presentation 

2.1. Disc structure 

Let us consider a planet on a circular orbit around a central star. 
In steady state, in the frame corotating and drifting along with 
the planet, the gaseous disc is split into 4 closed regions : 

1 . The Roche lobe of the planet, which is where the streamlines 
are closed around the planet. In the restricted three-body 
problem, it can be approximated by the Hill sphere, centred 
on the planet, of radius the Hill radius: r# = a p (q/3) 1 ^ 3 , 
which is where a p is the semi-major axis of the planet and 
q = Mp/M* is the ratio of the planet and stellar masses. 

2. The horseshoe region around the orbit of the planet, where 
the streamlines have a horseshoe shape. This region is delim- 
ited by the horseshoe separatrices. 

3. The inner disc between the star and the horseshoe inner sep- 
aratrix. 



4. The outer disc beyond the horseshoe outer separatrix. 

The Roche Lobe disappears for low -mass planets (q < 1CT 4 , but 
this depends on the aspect ratio, see lMasset et al. 2006, Fig. 11). 

The gas in the inner disc and in the outer disc exert a 
torque on the planet, the differential Lindblad torque. The gas 
present in the horseshoe region exer ts the so-called horse- 
shoe drag dWardl I199U iMassetl 1200 ll) . For low-mass planets 
in locally isothermal discs, the sum of these two effects is 
generally a neg ative torque and responsible for type I plane- 
tary migration (IWardl [1997). If the energy equation is taken 
into account, the torque due to the horseshoe drag is modi- 
fied in radiatively inefficient discs as well as in more realis- 
tic radiativ e discs. This may slow down or reverse the type I 
migrat i on dPaardekooper & M ellema 2006; Baruteau & Masset 
l2008at Kiev & Cridall2008l) . For intermediate mass planets, un- 
der some circumstances, th e horseshoe drag can also be respon- 
sible for type III migration dMasset & Papaloizoul 2003). 

The effect of the gas present in the Roche lobe of the planet 
has so far bee n neglected by most a u thors (with the notic eable 
exceptions of ID' Angelo et alj 120051 : iPeplihski et al.l 2008, see 
below), and never explicitly analysed. However, as this material 
is located very close to the planet, it could exert on it a strong 
force and torque. Consequently, it has to be carefully taken into 
account. We show further that it can influence the migration rate 
by almost a factor 2. This concerns planets of sufficient mass, 
which are able to create a Roche lobe in the disc and to accrete 
gas. 

Please note that the CPD does not necessarily fill the entire 
Hill sphere. There may be material inside the Hill sphere that 
is not orbiting around t he planet, but that is simply passing by. 
ID' Angelo etail d2005l) : iPeplinski et~aT] d2008l) have shown that 
the gas dynamics around a migrating planet is more complicated 
than described in Sect. 12.11 and that the Roche lobe may be mod- 
ified. This may be a source of confusion. Hereafter, the CPD al- 
ways refers to the material bound to the planet - that is: on a 
streamline that describes an orbit around the planet- while the 
Hill sphere refers to the material located within a distance r# to 
the planet. 

2.2. Gravitational interactions 

In Fig. Q] are displayed all the gravitational interactions phys- 
ically present in the problem. The gaseous disc has been split 
into the CPD, which is the gas orbiting the planet, bound to the 
planet, and the rest of the ProtoPlanetary Disc (PPD). Four sys- 
tems are interacting with each other, which makes a total of 12 
interactions. In addition, the PPD and the CPD, as extended sys- 
tems, exert on themselves a gravity force. 

These 14 gravitational actions can be classified in 4 cate- 
gories : 

- 5 direct terms (red solid arrows): The Planet— >PPD, 
Planet->CPD, Star->PPD, Star^CPD, Star->Planet are the 
most essential interactions of the system. They are responsi- 
ble for the disc and planet orbiting the star and for the per- 
turbations of the disc, including the mere existence of the 
CPD. 

- 3 indirect terms (red dashed arrows) : The (re)actions of the 
PPD, the CPD, and the planet on the star are generally con- 
sidered as indirect terms in the expression of the potential, if 
the frame is centred on the star. 

- 4 selfgravity terms (blue dotted arrows) : All the interac- 
tions between parts of the gaseous disc belong to what is 
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Fig. 1. Summary of the gravitational interactions in the problem 
of a planet, surrounded by a CPD, orbiting in a PPD, about a star. 



called the selfgravity of the disc (PPD^PPD, PPD^CPD, 
CPD^PPD, CPD^CPD). 
- 2 migration terms (green double arrows) : The action of the 
disc (PPD and CPD) on the planet are responsible for the 
planetary migration. 

If all the 14 interactions were taken into account in the sim- 
ulations, the computations would be self-consistent. The inter- 
actions of the two first categories are always taken into account 
and are clearly understood. But the selfgravity is generally not 
computed, because it is too expensive in terms of CPU-time. 
Consequently, the 2 migration terms may be inaccurate. 

The PPD< — >PPD effects, which is the influence of the gas 
selfgravity on the differ ential Lindblad torque, has been studied 
for low -mass planets bv lNelson & Benzl d2003l).|Pierens & Hurel 
(2005), and more recently by iBaruteau & Massell d2008blcf) . 
They confirmed that the type I migration is slightly accelerated 
by the disc selfgravity. This effect is not discussed here. On the 
contrary, we only focus on the effect of the CPD. 

In this paper, we plan to analyse the role and the behaviour 
of the CPD in the migration, which is represented by the arrows 
CPD^Planet, CPD->PPD, and PPD^CPD. As already said, the 
action Planet— >CPD should always be present and is not in ques- 
tion. In simulations with full gas selfgravity computed, all the 
blue dotted arrows of Fig.Q]are taken into account, including the 
CPD<->CPD arrow. However, the selfgravity of the CPD should 
only give a minor change of the CPD structure. It should there- 
fore be negligible for what concerns the planetary migration. In 
what follows, we do not discuss the role of the CPD<->CPD grav- 
itational interaction. 

2.3. Effect of the circumplanetary disc on the planet 

The presence of a CPD, gravitationally bound to the planet, can 
affect the migration by three means : 

1. CPD-> Planet: External perturbations (Star->CPD, 
PPD— >CPD) make the CPD not axisymmetric around 



the planet. Thus, the CPD could exert a torque on the planet ; 
we call this torque due to the CPD structure a direct torque. 

2. CPD->PPD->Planet : By perturbing the PPD, the CPD in- 
fluences the action of the latter on the planet. Basically, the 
CPD acts like the planet does, enhancing the amplitude of the 
perturbation. We define this as the perturbing mass problem. 

3. CPD— >Planet in migration: The planet and all the material 
that is bound to it migrate together. During the migration, 
the CPD is still linked to the planet. Consequently, the CPD 
is like a ball and chain for the planet. If the full gas self- 
gravity is computed, the CPD feels almost the same specific 
torque from the PPD than the planet ; then, the ball is pulled 
by the PPD, together with the planet, and the chain is not 
taut. If the full gas selfgravity is not computed (no blue ar- 
row in Fig. Q3, the CPD feels no torque from the PPD and 
tends to stay on a constant orbit. Then, the migrating planet 
has to pull the CPD. The chain is taut, which slows down 
the planetary migration. We define this as the inertial mass 
problem. 

The first and third cases are included in the green, double 
arrow CPD— >Planet in Fig.Q] The force acting on the planet in 
these two cases arises from gravitational interaction with the gas 
located in the CPD itself. 

The second case is an indirect interaction, represented in 
Fig.[TJby the two-arrow path: CPD— >PPD— >Planet. The action 
of the CPD on the PPD (blue dotted arrow CPD->PPD) may 
change the value of the green double arrow PPD— >Planet. 

2.4. Problem analysis and possible solutions 

In this paper, we evaluate either the specific torque felt by the 
planet T p or the migration rate d p (the dot denotes the derivative 
with respect to time). Both quantities are equivalent, as for a 
planet on a circular orbit, 2a p a p B = T p , with B the second Oort 
constant. 

It is well-known that a planet orbiting in a gaseous PPD per- 
turbs the latter gravitationally. In the linear regime, this perturba- 
tion leads to the appearance of a one-armed spiral wave, called 
the wake, in the inner and the outer discs. The amplitude of 
the perturbation is proportional to the planet mass M p . The per- 
turbed density distribution leads to a gravitational force on the 
planet. This force is proportional to the mass of the wake times 
the mass of the planet i.e. proportional to M p 2 . Thus, the planet 
feels from the disc a torque proportional to its mass squared. As 
a result, the specific torque applied to the planet is proportional 
to the planet mass M p . This corresponds to the PPDi^Planet ar- 
rows in Fig.Q] 

If the planet is surrounded by a CPD, things are slightly dif- 
ferent. Indeed, the CPD is gravitationally bound to the planet, 
and it moves with the planet. Thus, the planet of the above para- 
graph should be replaced by the planet and its CPD. The specific 
torque felt by this system is proportional to (M p + Mcpd), where 
Mcpd is the mass of the CPD. The presence of the CPD should 
increase the migration speed. 

More generally, the specific torque felt by the migrating body 
is proportional to the amplitude of the perturbation of the disc 
(proportional to the perturbing mass of the planet), times the 
mass of the body gravitationally interacting with it {gravitational 
mass of the planet), divided by the mass of the migrating body 
(inertial mass of the planet). If selfgravity is computed, these 
three masses are equal to (M p + Mcpd)- However, as we see 
below, the perturbing, gravitational, and inertial masses of the 
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planet differ if the gas selfgravity is not computed. This is the 
problem. 

In numerical simulations without gas selfgravity, the CPD 
does not perturb the PPD, neither does it feel a torque from it : 
in Fig.Q] the blue dotted arrows are suppressed. The perturbing 
mass is M p , and the force felt by the planet is proportional to 
M p 2 . But the CPD is linked to the planet. To follow the planet in 
its migration, the CPD has to exchange angular momentum with 
it. Thus, the relevant inertial mass of the planet is not M p but 
(M p + Mcpd)- The specific torque felt by the migrating system 
(the planet and its CPD) is proportional to M p 2 /(M p + M C pd)- 
Here, the presence of the CPD decreases the migration speed, 
which should not occur. 

Based on this analysis, several authors (e.g. 
iPierens & Nelson! 120081) have tried to cancel the force from 
the CPD on the planet by excluding the Hill sphere of the 
computation of the force of the disc on the planet. This aims 
at suppressing the green double arrow CPD— >Planet in Fig. Q] 
Thus, the inertial mass of the planet is M p again, and one is 
back to the first case ; the specific torque is proportional to M p , 
but not to (M p + Mcpd)- 

Others (e.g. iMasse J 120061: iPeplinski et al.l l2008h have sug- 
gested adding Mcpd to the p erturbing mass of the planet (de- 
noted M* by [Peplihski et al.h . to make the amplitude of the 
wake proportional to (M p + Mcpd)- In Fig. [1] this is equiva- 
lent to adding the blue, dotted CPD— >PPD arrow. This gives a 
torque proportional to (M p + Mcpd)M p and a specific torque 
proportional to (M g + M C pd)M p /(M p + M C pd) = M p as well. 
Pepli hski et al.l d2008l) find very little change in the migration 
rate by using M p or (M p + Mcpd) as the perturbing mass of the 
planet, which is quite surprising, as they have Mcpd ~ M p . 

By adding Mcpd to the planet mass felt by the disc (the 
perturbing mass) and excluding the CPD from the disc felt by 
the planet, one should recover a specific torque proportional to 
(M p + Mcpd)M p /M p = (M p + Mqpd)- One could also think of 
adding Mcpd to the perturbing mass of the planet and to the grav- 
itational mass of the planet. In that case, the torque is propor- 
tional to (M p + Mcpd) 2 - As the inertial mass is (M p + Mcpd), we 
find a specific torque in (M p + Mcpd), as required. 

Peplihski et al. (2008) have already addressed this issue and 
quite clearly described the problem. In the frame of type III mi- 
gration, they argue against any exclusion of a part of the disc 
in the torque calculation, but they propose another solution by 
applying the acceleration felt by the planet from the disc to 
the CPD (see their Eqs. (13) and (17)). This should mimic the 
blue, dotted arrow PPD— >CPD. In their paper, this acceleration 
is given by all the disc, CPD included. Then, if the CPD exerts, 
say, a positive torque on the planet for whatever reason, this pos- 
itive contribution will be applied to the CPD as well. This should 
not be, because it opens the possibility for the CPD to pull the 
planet indefinitely, without losing any angular momentum. For 
consistency, the acceleration applied to the CPD should be the 
one felt by the planet from the PPD, CPD excluded : the green, 
double arrow PPD— » Planet in our Fig.Q] 

The giant planets studied in this paper are not in the linear 
regime, so that the torque is not exactly proportional to M p 2 . 
However, this approximation allows the above reasoning. It en- 
lightens the role of the CPD in the migration process and sug- 
gests solutions. In Sects. [4] and [5] we numerically test these 
recipes and compare them with the results obtained with full 
selfgravity. 
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Fig. 2. The function fb(s) given by Eq. (01 for various values of 
b. 

3. Code and setting description 

To compare the migration rates obtained with different prescrip- 
tions to take the CPD into account, we performe d numerical sim- 
ulations, usingthe code FARGO dMasse t 2000a b). This is a pub- 
licly available^ 2D code, using polar coordinates (r, 6). We used 
the ADSG version, in which the energy equation and the gas 
selfgravity are implemented . The gas selfgravity is described in 
iBam teau & Massed (l2008d) . It can be turned on or off. 

The uniform kinematic viscosity is v = 10~ 5 a p 2 Q, p (Q, p be- 
ing the angular velocity of the planet and a p its semi major axis). 
The initial aspect ratio is ho = (H/r)o = 0.03. To take the disc 
thickness into account, a smoothing length e = 0.018 a p (60% 
of the initial disc height) is used in the expression of the planet 
potential $> p : 



<b P = - 



gm; 



(1) 



where s' = Vs 2 + e 2 with s the distance to the planet, and M* is 
the perturbing mass of the planet, equal to M p if not otherwise 
specified. 

The vector force F p j, exerted by the disc on the planet is 
computed using the following expression : 



<P,b 



Js=0 Jb=o 



2n GMpY. 



f b (s)s &6&s, 



(2) 



where s is the vector originating from the planet to the centre 
of the considered cell. The force F p j, depends of the parame- 
ter b through the term fb(s) : this is our filter used to exclude 
smoothl y the neigh bourhood of the planet, if necessary, also 
used in ICrida et all d2008l) . It is a smooth increasing function 
from at s = to 1 when s — » oo through 1 /2 when s - brn, 
drawn in Fig. [2] and given by: 



Ms) 



e X p|-10[— - 1|| + 1 



(3) 



For s < bra, M s ) ^ 1> while fi,(s) ~ 1 for s > brn, so that 
the gas within a distance br H of the planet is neglected while the 
gas outside is taken into account. For b — 0, fo = 1, so that all the 
disc is taken into account ; b — gives the standard expression 
of the force, with no exclusion of any part of the disc. 

The planet is not accreting. 

In the energy equation, the heat source is the viscous dissi- 
pation, and the cooling follows the following law, different from 
the public version of FARGO-ADSG : 



_ = -o- R T /kL , 
ot 

1 |http : //f argo . in2p3 . fr/| 



(4) 
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Table 1. Resolution of the simulations. 



time [orbits] 


resolution 

Nri,^ x Sectors Sr/r = 69 


designation 


from to 150 




115 x 393 


0.016 


initialisation 


from 150 to 250 




230 x 786 


8 x 10~ 3 


low 


from 250 to 325 




460 x 1572 


4 x 10~ 3 


middle 


from 325 to 400 




920x3144 


2 x 10~ 3 


high 


from 400 to 477 




1840x6288 


io- 3 


very high 



where e is the thermal energy of the gas per unit area, cr R is the 
Stefan-Boltzmann constant, T is the temperature, £ is the surface 
density of the gas, and k is the opacity, given by : 



k = <t r 



9S„ 2 v 



(5) 



This arbitrary constant opacity k ensures that, in the absence of 
perturbation, the initial state is in thermodynamic equilibrium 
with a density profile ~L(r) = Sor -1 ^ 2 . Our equation of state is 
then, with c s the sound speed, yu the mean molecular weight, R 
the universal gas constant, and y = 1 .4 the adiabatic index : 



(6) 



The computation of the energy equation provides more re- 
alistic discs. In addition, it reduces the problem of the accretion 
of the CPD. Indeed, with a locally isothermal EOS, no pressure 
effect stops the collapse of the CPD. Then, at high resolution, 
gas in the neighbourhood of the planet reaches extremely high 
densities, and convergence is hard to obtain. Since we plan to 
increase gradually the resolution of our simulations (see next 
paragraph), we need a viscous and a compressional heating and a 
radiative cooling for the CPD structure to remain within reason- 
able limits of density and temperature. In addition, preliminary 
tests show that with the gas's full selfgravity and locally isother- 
mal EOS, we could not reach convergence of the mass inside the 
Hill sphere. This motivated us to treat the gas thermodynamics. 
Our opacity is admittedly not realistic, and another prescriptions 
for k would make the CPD larger or smaller. Once again, we are 
not interested in the true structure of the CPD in this study, but 
we simply want a stationary CPD, of reasonable size, to check 
how it influences the migration. 

The radial grid extent is r e [Q.4a p ; 2.5a p ]. The regions 
[0Aa p ; 0.5a p ] and [2. 1 a p ; 2. 5a p ] are wave-damping regions, in 
which the density, the temperature, and the velocities are damped 
toward the initial values. The rings are geometrically spaced and 
the cells squared: 6r/r = 66 is constant. To save computing 
time, the resolution is initially low, and regularly doubled (see 
TableQ}- In the end, a resolution of A^ngs xN sect0K = 1840x6288 
and 6r/r — 66 = 10 -3 is reached. To our knowledge, this is the 
first time that the FARGO code is run at such a resolution. Most 
of the very-high-resolution simulations were run on a parallel 
cluster with 64 CPUs. The FARGO code uses a staggered mesh, 
with the density and temperature evaluated at the centre of the 
cells, and the velocities at the boundaries. When the resolution 
is doubled, each cell is divided into four. The new values of the 
fields were computed with a 2D linear interpolation on the old 
grid. 

Some remarks concerning the gas selfgravity are in order. 
The ADSG version of the FARGO hydro-code enables switch- 
ing on or off the gas selfgravity, the co mputation of which is de- 
scribed in Baruteau & Masset (l2008d) . In addition, it is possible 
to compute only the axisymmetric component of the gas gravity. 



This is the gravity field caused by the azimuthally averaged den- 
sity field of the gas, or by an axisymmetric disc with the same 
density profile as the actual gas disc. This component of the gas 
gravity field is responsible for increasing the angular velocity 
of the planet and of the gas itself (at least far enough from the 
disc inner edge). If one abruptly switches on the gas selfgrav- 
ity when restarting a simulation computed without selfgravity, 
a discontinuity in the velocities is introduced. To avoid this, all 
our simulations are computed with the axisymmetric component 
of the gas gravity taken into account. This has a negligible com- 
putational c ost, while the full se l fgravity is quite expensive. As 
discussed in Baruteau & Masset (2008c), accounting for the ax- 
isymmetric component of the disc gravity has only a marginal 
effect on the Lindblad torque with respect to the situation where 
the planet is held on a fixed circular orbit in a non selfgravi- 
tating disc. In the cases that we study here, the planet is mas- 
sive enough to shape a (deep or shallow) gap in the disc and to 
have a CPD. Consequently, the effect of selfgravity is different 
from in the above-mentioned studies. The shift of Lindblad reso- 
nances -which is essentially the only effect in type I migration - 
is taken over by the question of the gap shape, which gives the 
gas density at the position of each resonance, hence the torque. 
In addition, there are other issues, such as a modification of the 
perturbing, gravitational, and inertial mass of the planet, which 
we will address here. 



4. Type II migration case 

Following the procedure described in Sect. [3] we performed sim- 
ulations of a Jupiter mass planet, with gas initial surface den- 
sity: Eo(r) = 10~ 4 (r/a n )~^ 2 M*/a D 2 . In this case, the plane t 
opens a deep gap dLin & Papaloizou 1986a; ICrida et aT1 l2006). 
and the disc mass is low enough for type III m igration not 
to happen (Fig. 12 of Masset & Papaloizou] 2003b. Thus, the 
planet should migrate in type II migration dLin & Papaloizou! 
1986b). This choice could look inappropriate, because in stan- 
dard type II migration, the migration speed is independent of 
the planet mass. However, this stands only if the disc is massive 
enough to push the planet at the viscous accretion speed. This 
occurs when n > q/4, where u - Ana D 2 l.olM, and q - M p /M, 
dCrida & Morbidellill2007t ICrida! 1200 6). Here, q > 4//, the disc 
is not massive enough for the migration to occur at a viscous 
rate. Consequently the migration rate should decrease with the 
inertial mass of the planet. This case is timely for studying our 
problem. 



4. 1 . Inertial mass problem 



As explained in Sect. I2.3I this problem only concerns migrating 
planets. To study this effect, we let the planet evolve freely in 
the disc. To test and quantify the influence of the CPD on the 
migration speed, a region around the planet is excluded from the 
computation of the force of the disc on the planet. To perform 
this, the parameter b in Eq. (O is taken strictly positive, and sev- 
eral values are tested. 



4.1.1. Middle resolution 

First, the planet is released after 300 orbits on a fixed circular 
orbit ; the resolution of the grid is 6r — 4 x 10~ 3 r = rn/17.3 (see 
Table [1), which is finer than the typical kind of resolution used 
in the literature. The migration is followed for 25 orbits. Various 
values of b from to 1 are used. The results are displayed in 
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Fig. 3. Migration path of the Jupiter mass planet after its release 
at t = 300 orbits, with 6r/r - 69 - 0.004. The bold curve cor- 
responds to a simulation where the full gas selfgravity is com- 
puted. The other ones come from simulations with exclusion of 
a part of the Hill Sphere, for various values of b. The first five 
orbits are enhanced in the inset at the bottom left. 



at the moment where the planet is released, computed using the 
following expression : 

M H (b)= E (I -Ms)) sd9ds. (7) 

Js=o Je=o 

If, instead of the smooth ft, function, a Heaviside function was 
used, Mu(b) would have the more intuitive following expres- 
sion : 

pbr H pin 

Mn(b) ~ E sd9ds . 

Js=o Je=o 

If E was constant, Mu(b) should increase as b 2 . This is approx- 
imately the case only for b < 0.15, which shows that the den- 
sity decreases strongly with s for s > 0.15r#. In fact, 90% of 
the mass is inside 0.55r//, which favours the idea that the CPD 
should be considered as only the inner half of the Hill sphere. 

The total mass in the Hill sphere is 0.61% of the planet mass, 
so that the inertial effect should be very small. In addition, the 
variation of Mn{b) for b > 0.5 is only a few percent, so the 20% 
difference in the migration rates between the cases b - 0.5, 0.6, 
0.8, and 1 .0 cannot be explained by a variation in the mass of the 
excluded region. 



0.007 




Fig. 4. Mass distribution inside the Hill sphere of the Jupiter 
mass planet given by Eq. @, in planet mass units. The mea- 
sure is made after 300 orbits on a fixed circular orbit, with a 
resolution of 6r/r = 0.004. 



Fig. [3] The curves for b < 0.3 almost overlap ; therefore, only 
the curve for b — is drawn as a thin red solid line with + 
symbols. For b = 0.4, the planet follows the blue starred line. 
It migrates significantly more slowly, and in 25 orbits, its semi- 
major axis has decreased 8.5% less. For b > 0.4, the migration 
speed increases with b, and for b = 1 (exclusion of all the Hill 
sphere), the semi-major axis of the planet decreases in 25 orbits 
nearly 12% more than in the case b — 0. The migration ratio after 
10 orbits between b = 0.4 and b = 1 is almost 2. This shows that 
the Hill sphere can have a strong influence on the migration rate 
in numerical simulations. 

The bold solid red curve in Fig. [3] is obtained with a simu- 
lation where the full selfgravity of the gas has been taken into 
account, while only the axisymmetric component is used in the 
others. Of course, in that case, no exclusion of any part of the 
disc is done, b — 0. This bold curve should be seen as the refer- 
ence. 

To better understand what happens in the Hill sphere, Fig. [4] 
displays the mass distribution in the Hill sphere of the planet, 



4.1 .2. High resolution 

In a second experiment, the planet is not released at t — 300, but 
kept on a fixed circular orbit until 375 orbits, the resolution being 
doubled at time 325 (see Table[TJ. The planet is then released at 
375 orbits, with a resolution of 0.002 (35 cells in a Hill radius, 
1000 cells in the surface of a circle of radius 0.51r#). The results 
are displayed in Fig. [5] The curves with symbols were performed 
using fb as filter, for various values of b. 

The average migration rate is a bit slower than in the previ- 
ous case. The green dashed curves with x symbols (case b — 1) 
in Figs.[3]and|5]are rather linear, providing a constant migration 
rate, only given by the torque from the PPD outside of the Hill 
sphere. This torque should be converged for resolutions of the or- 
der o f rn/10 (here: 0.0046 ), as widely assumed in the literature 
(e.g. iD'An gelo et af1 l2005l) . However, the gap profile is slowly 
evolving with time, so that the torque exerted on the planet at 
t = 300 or 375 orbits differs, independently on the resolution. 
This time delay is the reason why the green dashed curves with 
x symbols slightly differ in Figs.[3]and[5] 

For the other curves, with narrower filters fi, or full gas self- 
gravity, their various evolutions between the medium-resolution 
and high-resolution cases suggest that convergence in the Hill 
sphere had not been reached at medium resolution. In particular, 
the different variations of the cases b = 0.4 and b = 0.6 between 
Figs. [3]and[5]can hardly be explained by a simple smooth time 
variation of the disc structure. However, we recall that a conver- 
gence study of the CPD structure is certainly interesting but is 
not within the scope of this paper, focused as it is on the role of 
the CPD in migration in standard numerical simulations. 

With 6r/r = 0.002, all the curves are almost linear, denoting 
that the CPD structure is better resolved. Then, the role of the 
filter function can be analysed better. Here, the trend is clear. 
The larger b, the lower the inertial mass of the planet, and the 
faster the migration, as expected if excluding a part of the Hill 
sphere only had an effect on the inertial mass of the planet. 

In addition, other shapes of filter are tested. The black dots 
correspond to a case where fi,(s) has been replaced in Eq. (|2]i by 
a Gaussian filter : 

g(s) = 1 - exp(-(s/r H ) 2 ) . (8) 
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Migration path of the Jupiter-mass planet after its release 
375 orbits, with 6rjr = 69 — 0.002. Same colour code as 
The five first orbits are enhanced in the inset at the bottom 



is probably not due to the decrease in the inertial mass of the 
planet. As in the middle resolution case, we find that by exclud- 
ing all the Hill sphere, one misses a part of the torque that this is 
not an artifact from the absence of selfgravity. 



4.2. Perturbing and gravitational mass problem 

The addition of the CPD mass to the perturbing mass of the 
planet can be easily implemented in the code. One should sim- 
ply add the mass of the CPD to the planet perturbing mass in 
computating the gravity potential of the planet through Eq. (HJ. 
However, it is not easy to measure the mass of the CPD Mcpd, 
because this requires analysing the disc structure. We assume 
that the CPD is a disc centred on the planet of radius r « 0.6r#. 
Therefore, we take Mcpd = Mu(b = 0.6). Figures [4] and [6] have 
shown that the obtained value of Mcpd does not vary by more 
than 10% for 0.5 < b < 1. The CPD itself should not feel its 
own mass as added to the planet. This would be inconsistent. 
Therefore, the expression of the gravitational potential field of 
the planet is now given by Eq. (HJ, with 
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Fig. 6. Mass distribution inside the Hill sphere of the Jupiter 
mass planet given by Eq. (|7), after 375 orbits, with a resolution 
of 6r/r = 0.002. 

It gives almost the same result than fo g (light blue dot-dashed 
curve with full squares). The green dashed line without x sym- 
bols corresponds to a Heaviside filter h(s) = if s < rn, h(s) = 1 
if s > rn- It gives almost the same result as f\ (green dashed line 
with x symbols). 

For small b, the difference between b = 0.6, b = 0.4, b = 0, 
and the full selfgravity case is rather small, although no match 
is observed. The effect of gas selfgravity on type II migration is 
little and has nothing to do with its effect on type I migration, 
which could be expected as the torque felt by the planet in this 
case is equal to the viscous torque from the disc and not the 
differential Lindblad torque (which is affected by the gas full 
selfgravity). Higher values of b, or the use of g{s) or h(s) show a 
significant acceleration of the migration with respect to the full 
selfgravity case. 

Figure [6] shows the mass distribution inside the Hill sphere 
of the planet at the moment it is released, after 375 orbits. The 
distribution is a bit more compact than in previous case : 90% 
of the mass is now included in 0.48r//. The total mass is higher 
(0.01 M p ). This is due to the doubling of the resolution. 

The value of Mu(b) does not vary significantly for b > 0.5. 
Consequently, the inertial mass of the planet should not change 
much whether one half of the Hill sphere is excluded or the 
full Hill sphere. The acceleration of the migration for b > 0.6 



M; = M p + Mcpd X JiCO , (9) 

so that inside brn, the potential is basically unchanged. Outside, 
the effective perturbing mass of the planet is M p + Mcpd. 

Ideally, for s < brn, one should use Mh(s/Vh) as defined by 
Eq. (0 instead of MqpdM s )- B ut this is computationally more 
expensive as it requires calculating Mu(s/rn) for every cell of 
the grid within a distance brn to the planet. Alternatively, one 
could also use a function other than fy, which would fit the 
curves displayed in Figs.|4]and|6]better. For consistency, we have 
kept fb. This should have only a minor influence on the CPD, and 
even less on the global simulation. However, it is important that a 
smooth enough function is used : with a Heaviside-like function, 
(S)p(s) could not be monotonic. In the present case, / ' 6 (s) < 4.2, 
while M C?D /M p * 0.01, so that the change in <b' p {s) at s = 0.6r H 
does not exceed 4%. 

With this prescription, we restarted our simulation from t - 
350 orbits, with the planet on a fixed circular orbit, and a high 
resolution of 6r/r — 60 — 2 x 10~ 3 . The mass of the CPD is 
added smoothly over 3 orbits to the planet mass to avoid a dis- 
continuity : in Eq. ©, Mcpd is replaced by Mcpd * g(t) with 
g(f) = sin 2 ((f - 350)tt/6) for t < 353 and g = 1 for t > 353. 

At 375 orbits, the planet is released and free to migrate, like 
in Sect. 13X21 At t his tim e, the CPD mass is 1.02 x 10~ 5 M, = 
\%M p , like in Sect. 14.1. 21 The result in terms of migration is dis- 
played in Fig.Qas the curve labelled "large M* ; b=0". The thin 
red solid curve labelled "M* = M p ; b — 0" is the most stan- 
dard simulation, where nothing is done (no full selfgravity, no 
exclusion, no modification of the masses), also present in Fig. [5] 
Surprisingly, the migration is slowed down in the first three or- 
bits when the perturbing mass of the planet is increased. The 
green curve with x symbols corresponds to a case where M* is 
given by Eq. ©, and in addition b = 0.6 in Eq. The black 
curve with big dots corresponds to a simulation where the per- 
turbing mass is given by Eq. (0 and the gravitational mass of 
the planet is M p + Mcpd- In the first five orbits, adding Mcpd to 
the gravitational mass of the planet or excluding 60% of the Hill 
sphere both lead to a similar acceleration of the migration : the 
green and black curves overlap at t — 380, below the blue curve. 
However, in the longer term (20 orbits), these two recipes lead 
to different migration paths. 



8 



Crida et al. : Dynamical role of the circum-planetary disc 



0.9995 



0.999 



0.9985 



0.998 




375 



380 



385 390 
time [orbit] 



395 



400 



Fig. 7. Migration path of the Jupiter mass planet after its release 
at t — 375 orbits. The label "large M*" denotes simulations for 
which the planet perturbing mass is given by Eq. (0. The label 
"large M^ 8 "™" means that the gravitational mass of the planet is 
M p + Mcpd- The bold curve is the reference case, with full gas 
selfgravity (and M* = M p ). 



as discussed in Sect. 12. 4| : instead of F p> o/M p , F Pi o,s/M p is ap- 
plied to the CPD, but the planet still feels F p q. During the first 
15 orbits, this has only a marginal influence on the migration 
and the two curves overlap. But after t = 390, the paths dif- 
fer strongly. This proofs that our reasoning in Sect. 12. 41 was not 
vain : it is important to decide whether the acceleration applied 
to the CPD should include the acceleration of the CPD on the 
planet or not. The path in the second case is closer to the bold 
red reference full selfgravity case. We recommend that the ac- 
celeration of the CPD on the planet should not be applied to the 
CPD itself. 

The green dashed curve with x symbols is taken from Fig. [7]: 
M* p is given by Eq. © in Eq. ®, and b = 0.6 in Eq. ©. The 
green and orange curves are very close but not identical. The 
acceleration of the CPD or its exclusion of the force exerted by 
the disc on the planet both give the planet an inertial mass of M p , 
but in a different way. 

Accelerating the CPD is appealing, because this method 
gives the planet an inertial mass of M p , without having to ex- 
clude a part of the disc. Moreover, it can be used with the addi- 
tion of Mcpd to the perturbing mass of the planet. However, it 
does not give the same migration path as the gas full selfgravity 
either. 
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Fig. 8. Migration path of the Jupiter mass planet after its release 
at t = 375 orbits. The perturbing mass of the planet is given by 
Eq. © in every case, except the full selfgravity case (bold solid 
red line). The acceleration applied to the CPD varies. 



4.3. Acceleration of the CPD 



Last, we tested the recipe proposed by Pephris ki et alj (|2008). 
The perturbing mass of the planet is given by Eq. (0. The force 
felt by the planet is F p q, given by Eq. (0. The acceleration 
(or the specific force) felt by the planet from the disc is then 
FpA)/M p . This acceleration is then applied to the material in- 
side the Hill sphere, after multiplication by 1 - fo.s(s), so the 
additional acceleration is smoothly vanishing toward zero when 
s > 0.5 rn- The material inside 0.5 r# feels the same acceleration 
exactly as the planet and should naturally follow the same orbit 
and the same migration path. This should set the planet free of 
the "ball and chain" effect of the CPD. 

The resulting migration path is shown in Fig. [8] as the red 
curve with + symbols. The orange curve with open triangles is 
computed with the same algorithm, except that the acceleration 
felt by the CPD is the one exerted on the planet by the PPD only, 



4.4. Summary of the type II migration case 

We observed the migration rate of a Jupiter mass planet, using 
various numerical recipe to take the CPD into account. The mass 
of the CPD is only 1 % of the planet mass here, but we find that 
the way the CPD is considered can change the migration signif- 
icantly. In particular at lower resolution, excluding a part of the 
Hill sphere can lead to variations of almost a factor two. 

On the other hand, changing the perturbing or gravitational 
masses of the planet or accelerating the inner half of the Hill 
sphere leads to a moderate change in the migration speed. The 
resulting migration path is compatible with the full selfgravity 
case. Unfortunately, none of the methods leads to the same mi- 
gration path as when the full selfgravity of the gas is computed. 

This suggests that the mass of the CPD was in this case too 
low to strongly influence the migration, but that the outer layers 
of the Hill sphere are able to exert a strong torque on the planet. 
Even if the Hill sphere contains a small amount of gas, it should 
therefore be considered with care. This is be discussed in more 
detail in Sect. [6] 



5. Type III migration case 

The second case that we studied is a Saturn mass planet with 
£o(r) = 10~ 3 r~'^ 2 M»/flp 2 . This density is one order of magni- 
tude higher than in previous case, so that type III migration is 
expected for this planet. Indeed, except for what concerns the 
slope of the density profile, this c orresponds to the case stud- 
ied by Masset & Papaloizou ( 2003), who found type III migra- 



tion, and iD'Angelo et al 



"(120051) . who claim that type III mi- 
gration does not happen at high resolution. Another difference 
with these previous works is that our equation of state is not lo- 
cally is othermal, which could significantly change the corotation 
torque dBaruteau & Massetll2008at iKlev & C rida 2001). 
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Fig. 9. Migration path of the Saturn mass planet after its release 
at t — 300 orbits. Bold curve : simulation where the gas selfgrav- 
ity is computed. Curves with symbols : simulations with exclu- 
sion of a part of the Hill Sphere, using the filter fa, for various 
values of b. Curves without symbols : exclusion of a part of the 
Hill Sphere with Heaviside filters. Big dots : exclusion of the Hill 
sphere with a Gaussian filter. 



5.1. Inertial mass problem 

The same experiment as in the Jupiter mass planet case was 
performed. The setup was identical, with only the M p and So 
changed. 

5.1 .1 . Middle resolution 

For a restart at 300 orbits and resolution of 4 xl0~ 3 a p = r^jWA, 
the results are displayed in Fig. [9] The bold, red curve corre- 
sponds to the case where the full selfgravity of the gas has been 
computed. In the other cases, only the axisymmetric component 
of the selfgravity is computed. The curves with symbols come 
from simulations with exclusion of a part of the Hill Sphere, 
using the filter ft,, for various values of b (b — means no ex- 
clusion). For b < 0.3, the curves overlap and only the case b = 
is drawn. Then, for increasing b, the migration speed decreases. 

The black dots correspond to a case where fb(s) has been 
replaced in Eq. (0 by the Gaussian filter g(s) given by Eq. (O. 
It gives almost the same result than /0.9 (grey blue short dashed 
curve with full triangles). 

The green dashed line without symbols corresponds to ex- 
clusion of the Hill sphere with a Heaviside filter h{s) — if 
s < rn, h(s) = 1 if s > Th- It gives the most divergent result with 
the full selfgravity case. The pink dotted line without symbols 
corresponds to exclusion of the inner 60% of the Hill sphere 
with a Heaviside filter ho&(s) = if s < 0.6 rg, ho^{s) — 1 if 
s > 0.6 r/j. It gives a result similar to fo.j (orange triple dashed 
line with open triangles). Not only does the size of the excluded 
region matters but also the shape of the filter. 

First of all, we find that type I II runaway migration h appen s 
in that case, in agreement with Masset & Papaloizou (2003), 
even with full selfgravity (bold red line ). Admittedly however , 
we have not reached the resolution of ID' Angelo et al.l (2005) 
here (see next section). 

Second, when varying b, the migration rate changes dramat- 
ically, and for b ~ 1, the runaway migration process seems to 
be broken. This again shows that excluding a part of the Hill 



sphere or not can have strong consequences on the migration. At 
t = 325 orbits, the semi-major axis of the planet is a — 0.77 in 
the case b — and 0.93 in the case b = 1 . The ratio of the semi 
major axis variations is 3.33. 

Third, the observed behaviour is opposite the expected one 
if one considers only the inertial mass effect : here the migration 
speed decreases with b. The total mass of gas included in the Hill 
Sphere at t = 300 orbits is 2 x 10~ 5 M, = M p /16. The impact 
on the migration speed should thus not be more than 7%. This 
supports the idea that the observed difference in the migration 
rate when varying b is not due to a change in the inertial mass 
of the planet, but to the direct torque repartition inside the Hill 
sphere. The mass and torque distribution inside the Hill sphere 
will be discussed in Sect. [6] 

5.1.2. Very high resolution 

This simulation is run until t = 475 orbits, with the planet on a 
fixed circular orbit, and increasing resolution following TableQ] 
The axisymmetric component of the gas gravity is always taken 
into account. In the end, a resolution of 10~ 3 is reached, so that 
the length of a Hill radius is covered by 48.2 cells. This is al- 
most e xactly the settings of the run 2D5Gb of ID' Angelo etal] 
d2005h . for which convergence in resolution was reached. The 
planet is then released and let free to migrate, excluding a bigger 
or smaller part of the Hill sphere. The resulting migration paths 
are displayed in Fig.[10]for various filters. 

The migration speed is much slower than at m i ddle resolu- 
tion (Pig. |9j, as already claimed by D' Ange lo et al.l (12005) (note 
the different y-axis scale), but the migration rate exponentially 
increases, in a runaway typical of type III migration. In fact, the 
difference can arise from two independent effects : (i) as already 
pointed out for the type II migration case, the disc density profile 
has evolved with time between t = 300 and 475 orbits, indepen- 
dently of the resolution (see Fig. ITTb . Therefore, the torque felt 
by the planet at the release date differs. If it is smaller initially (as 
it appears to be from Fig.fTTTi. the torque remains smaller during 
the runaway, (ii) The e-folding time of type III migration scales 
as CMD/{CMD - M p ), where CMD is the coorbital mass deficit, 
and M„ should be understood a s the inertial mass of the planet 
(see iMasset & Papaloizou 2003j). Here, the total mass of gas in 
the coorbital region in the unperturbed disc is about 6 x 10~ 4 M», 
and the value of CMD is of the order of M p . Therefore, a little 
increase in the inertial mass of the planet (due to the increase in 
Mcpd) results in a significant increase in the e-folding time of 
the runaway migration. We think that this explains most of the 
difference between Figs. [9] and [10] The latter effect may also be 
responsib le for the vanishing o f type III migration at high reso- 
lution in lD' Angelo et al.l(l2005h (Fig. 8) because the mass inside 
the Hill sphere increases with resolution up to a few M p , accord- 
ing to their Fig. 9. This shows once again that the Hill sphere 
should be considered with care. In the type III migration regime, 
convergence in the Hill sphere structure had better be reached. 
Here, however, we are not looking for the real migration rate, but 
for the influence of the gas in the Hill sphere on the migration 
rate. This experiment is suitable for this purpose. 

The major change observed does not stem for a strong mod- 
ification of the disc structure : the density profile is drawn in 
Fig. [TT]at the two times where the planet is released. The gap 
width and depth are similar. The CPD is almost twice as mas- 
sive at very high resolution than at moderate resolution, though, 
reaching 1% of the planet mass. This promotes the idea that the 
material in the Hill sphere or the CPD plays a major role in the 
type III migration rate. 
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Fig. 10. Migration path of the Saturn-mass planet after its release 
at t — 475 orbits. Bold red curve : the full gas selfgravity is com- 
puted. Other curves : only the axisymmetric part of the gas self- 
gravity is computed. Thin solid red curve with + symbols : the 
whole disc is taken into account. The other curves come from 
simulations with exclusion of a part of the Hill sphere, with the 
filter fi, for various values of b, for the Gaussian filter g, or for 
the Heaviside filter h. 
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Fig. 11. Density profiles at t = 300 orbits (Sr/r = 69 = 0.004) 
and at t — 475 orbits (5rjr — 59 — 0.001), with a Saturn mass 
planet on circular orbit at r = I. 



In the previous section, the results showed that a large part 
of the torque responsible for the runaway migration originated 
in the Hill sphere. Therefore, the lower migration rate observed 
in a very-high-resolution case leads to the conclusion that the 
torque exerted on the planet by the material inside the Hill sphere 
decreases with resolution. 

However, for the topic of this paper, we see that even at very 
high resolution, excluding more than 60% of the Hill sphere has 
a significant influence on the migration speed. The larger the 
excluded area, the more the migration diverges with respect to 
the full selfgravity case. The Heaviside filter to exclude all what 
is beyond r# of the planet appears to be the least appropriate. 
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Fig. 12. Migration path of the Saturn mass planet after its re- 
lease at t — 300 orbits, for different perturbing masses. The la- 
bel "large M*" denotes simulations for which M* is given by 
Eq. © since t = 290. The bold curve labelled "Self Gravity" is 
computed with full gas selfgravity. 



5.2. Perturbing and gravitational mass problem, acceleration 
of the CPD 

5.2.1. Middle resolution 

With the same prescription as in Sect. I4.2I we restarted our sim- 
ulation from t = 250 orbits, with the planet on a fixed circular 
orbit and a resolution of 0.004. At this time, M C pd = M H (0.6) 
has reached 10~ 5 M» = 0.035 M p . It is added smoothly over 40 
orbits to Mp, to avoid a discontinuity. 

Then, we let the planet evolve freely under the influence of 
the disc at time t = 300 orbits. The migration path is shown in 
Fig.[T2]as the red thin solid curve with triangles (labelled "large 
M* ; b=0"), compared to the migration rate without taking the 
CPD mass into account (thin red solid curve, labelled "b=0"), 
and to the selfgravity case (bold solid curve). By excluding in 
addition the inner 60% of the Hill sphere (i.e. taking b = 0.6 in 
Eq. one gets the light blue, dot-dashed line with full squares 
(labelled "large M* ; b=0.6"). It can be compared to the similar 
case, with M* = M p , which is the pink, dotted line with open 
squares (labelled "b=0.6"). The green dashed curve with x sym- 
bols (labelled "large M* and M^ lav ") corresponds to the last case 



described in Sect. I2.4I : the perturbing mass is given by Eq. (0, 
and gravitational mass of the planet is M p + Mcpd- 

As expected, the migration is accelerated by the augmenta- 
tion of the perturbing mass of the planet. Adding Mcpd to the 
perturbing mass of the planet leads to a rise in the migration 
rate of 4% in both cases b — and b = 0.6. For the 15 first 
orbits, the curves a p (t) remarkably overlap if one replaces t by 
300 + (f - 300)/ 1.04 in the cases with M* = M p . Another 4% 
increase in the migration speed is observed when Mcpd is also 
added to the gravitational mass of the planet. This is consistent 
with the mass of the CPD. 

The "Self Gravity" simulation was launched at 300 orbits 
from the same disc as the other simulations here : before t = 300 
only the axisymmetric component of the selfgravity is com- 
puted, and the disc is perturbed by a planetary potential given by 
Eq. (Q} with M* given by Eq. The migration path is almost 
indistinguishable from the one in the previous section, restarting 
from the classical potential perturbed disc. 
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Fig. 13. Migration path of the Saturn mass planet after its release 
at t — 300 orbits. The label "large M*" denotes simulations for 
which M* is given by Eq. (J9]l since t = 290. The bold curve 
labelled ' Self Gravity" is computed with full gas selfgravity. 



Unfortunately, as the full selfgravity case gives a slower mi- 
gration than the standard simulation, adding the mass of the CPD 
to the perturbing and gravitational masses of the planet gives an 
even more divergent result. This means that in this case, the self- 
gravity plays another role as simply modifying the perturbing, 
gravitational, and inertial masses of the planet. 

As in Sect.gJ] the algorithm of Pepliri ski et al.l(l2008l) is also 
tested and compared to the other ones. The migration obtained 
with it, with also M* given by Eq. (0, is shown in Fig. Q~2] as 
a blue curve with stars, labelled "large M* ; acpr>=F p fl/M p ". It 
follows the path of the case with simply M* given by Eq. ||9). 

5.2.2. Very high resolution 

The same experiments are computed at very high resolution from 
t — 475 orbits. At this time, the mass distribution in the Hill 
sphere of the planet is shown in Fig. [16] and Mqpd = Af#(0.6) = 
1.933xl0~ 5 M„, = 0.068M P . For the cases with M* = M p +M cm , 
the simulation is restarted with this perturbing mass at t — 450 
orbits, with the planet on a circular orbit until 475 orbits, where 
the planet is released. As expected, the addition of Mcpd to M* 
accelerates the migration. This can be seen by comparing the 
thin solid red line to the thin solid red line with open triangles 
in Fig. [13] These two cases have b — but a different perturbing 
mass. During the 5 first orbits, the migration paths overlap ex- 
actly if one replaces t with 475 + (f- 475)/ 1.068 in the M* = M p 
case. In the longer term, however, the acceleration of the migra- 
tion is closer to 10% than 6.8%. 

With b = 0.6, the influence of the perturbing mass can be 
seen by comparing in Fig. Q~3] the pink dotted line with open 
squares with the light blue dot-dashed line with full squares. The 
migration rates compare exactly like in the b — case : 6.8% dif- 
ference during the 5 first orbits, 10% in the longer term. 

An acceleration of F p o.5/M p is then applied to the CPD 
(green dashed curve with x symbols, labelled "«cpd = 
FpO.s/Mp"). This leads to an overall acceleration of the migra- 
tion of about 7% with respect to the standard, b = case. This 
could be expected, but is very different than excluding the CPD 
from the computation of the force of the disc on the planet (pink 
dotted curve with open squares). However, both methods aim at 



delivering the planet from the inertia of the CPD. This shows 
that they are not equivalent. 

If in addition to the CPD acceleration, the perturbing mass is 
M* = M p + Mqpd (blue dashed curve with stars), the migration 
is accelerated again. The speed-up is exactly 6.8% with respect 
to the previous case (acceleration of the CPD with M* = M p ) 
during the first 5 orbits and about 8% in the longer term. Last, in 
the case where M* = M p + Mcpd and b = 0, the acceleration of 
the CPD leads to a speed-up in the migration of about 5%. 

5.3. Summary of the type III migration case 

The outer layers of the Hill sphere (s > 0.6rn) play a big role in 
the type III migration regime. Excluding more than 60% of the 
Hill sphere leads to a decrease in the migration rate and even to a 
suppression of the runaway phenomenon. In the full selfgravity 
simulations, the runaway is present, even at very high resolu- 
tion, although the migration rate is decreased. Adding Mcpd to 
the perturbing and/or gravitational masses of the planet clearly 
leads to an acceleration of the migration at the expected rate. 
Accelerating the CPD appears to have less of an effect at middle 
resolution, while it plays a significant role at very high resolu- 
tion. 

At middle resolution, no recipe gives a perfect match with 
the case where the full selfgravity of the gas is computed. We 
should also conclude that the gas selfgravity plays a more com- 
plex role in the type III migration than simply accelerating the 
CPD or letting the CPD perturb the PPD. 

At very high resolution, excluding a part of the Hill sphere 
also has a strong influence on the migration path. With b = 
0.6, the agreement with the full selfgravity case is very good. 
However, this cannot be explained by the analysis of Sect. 12.41 
according to which taking b = 0.6 should accelerate the migra- 
tion while it is slowed down here. This agreement is therefore 
most likely a coincidence. In any case, the type III migration 
rate is strongly reduced with increased resolution. 



6. Direct torque and Hill sphere structure 

6.1. Direct torque from the CPD 

A short analysis of the direct torque from the CPD, defined as 
the material bound to the planet, shows that it should be neg- 
ligible. It has been observed in numerical simulations with grid 
refinement that two spiral arms form in the CPD, because of tidal 
effects from the star. Thus, the CPD should not be considered as 
an axisymmetric structure, but it looks centro-symmetric any- 
way. At lower resolution, the spiral arms disappear; however, 
the axisymmetry of the structure is not perfect. Only the pertur- 
bation from the axisymmerty can be responsible from a torque 
on the planet. 

It would be absurd that the CPD alone directly exerts a torque 
on the planet and makes it migrate on the long term, because the 
CPD is linked to the planet ; this torque is an internal interaction 
in the closed {CPD+Planet} system. This system cannot change 
its orbit by itself. The only possible angular momentum trans- 
fer is between the spin angular momentum of the CPD about the 
planet and the orbital angular momentum of the {CPD+Planet} 
system about the star, similar to the slow down of the Earth ro- 
tation together with increase of the Moon-Earth distance. If one 
assumes for simplicity that the CPD is a constant density disc of 
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Fig. 14. Mass Mn(b) (solid line, Eq. (01) and Torque (dashed 
line) distribution in the Hill sphere of the Saturn mass planet at 
time t = 300 orbits. 



radius rcpo in Keplerian rotation about the planet, its spin angu- 
lar momentum is 

2nsZ ^/GM p ~s ds = -M CPD V^Mp V^d~ ■ (10) 
o J 

The orbital angular momentum of the planet is J p = 
M p y/GM* -\fa~p~. The radius of the CPD is generally of the or- 
der of 0Ar H . Thus JqpdIJ,, ~ 0.5(M CPD /M / ,)^ 2/3 . Even in the 
case of a CPD as massive as the planet and of a 3 Jupiter mass 
planet, the ratio would be inferior to 1%. In the cases studied in 
Sects. |4] and [5] one finds Jcpo/Jp ~ 5 x 10~ 4 . Thus, the spin 
angular momentum of the CPD is negligible with respect to the 
orbital angular momentum of the planet, and consequently the 
tidal effects of the star (or the PPD) on the latter can only have a 
negligible effect on the planetary migration. 



6.2. Mass and torque distribution inside the Hill sphere 

We consider the situation after 300 orbits in the case of a Saturn 
mass planet in a massive disc for the middle resolution. FigurefBl 
displays the mass inside the Hill sphere Mn(b). It increases sig- 
nificantly with b even for b > 0.4, because the gap opened by 
the planet is not empty (see Fig. [TTJ ; the background density 
is responsible for the increase in Mu with b. However, the az- 
imuthally averaged profiles shown in Fig.QT]show that a CPD is 
present. 

The torque exerted on the planet by the region within s < 
brn, is also plotted in Fig.[14]as a function of b. More precisely, 
it is the torque of the force computed with Eq. ©, with (1 -fb(sj) 
instead of fb(s) ; or equivalently the torque of the force : F p o - 
Fpj,. It appears that for a planet on a fixed circular orbit, the 
region inside r#/2 exerts a negligible torque, while the region 
between r#/2 and r# has a non negligible effect on the planet: 
the direct torque is not small. 

Figure [TBJ shows the torque exerted on the planet by vari- 
ous regions of the disc, as a function of time, during the type III 
migration between 300 and 325 orbits. The measures are done 
in the simulation with full gas selfgravity computed. The total 
torque, exerted by all the disc, is increasingly negative, which 
is characteristic of runaway inward migration (solid line). The 
other curves are the torques due to spheres centred on the planet 
inside the Hill sphere, like in Fig. [14] In the bottom panel, the 
torques due to fractions of the Hill sphere are plotted in terms of 
percentage of the total torque. The torque exerted by the material 
inside the Hill sphere is positive before the planet is released, but 
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Fig. 15. Torque exerted on the planet as a function of time, dur- 
ing its migration, by all the disc (bottom, solid line), the Hill 
sphere (dashed line b - 1.0), the inner 80% of the Hill sphere 
(b = 0.8), 60% (b = 0.6), and 40% (b = 0.4). Bottom panel : the 
same torques as a percentage of the total torque. 
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Fig. 16. Mass M H (b) (solid line, Eq. ©) and torque (dashed line) 
distribution in the Hill sphere of the Saturn mass planet at time 
t = 475 orbits. 



then turns negative. The inner 40% of the Hill sphere plays a neg- 
ligible role for a migrating planet too. But the entire Hill sphere 
is responsible for about 40% of the total torque. This confirms 
that the region between s = r#/2 and s — rn plays a major role 
in the process of type III migration, as Fig. [9] suggests. 

Figures [16] and [T7] are the same as Figs. [T4l and [TBI but from 
time 475 orbits on, with the very high resolution. The torques 
as a function of time in Fig. [17] come from the simulation with 
only the axisymmetric component of the selfgravity taken into 
account, and no exclusion of any part of the Hill sphere. The 
two figures look very similar to the previous ones, in particular 
the total torque is smaller but is also increasingly negative. The 
other previous remarks also apply at very high resolution. 

6.3. Size of the circum-planetary disc 

The calculation results and the analysis of the mass and torque 
distributions inside the Hill sphere suggest that the region within 
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Fig. 17. Torque exerted on the planet as a function of time, dur- 
ing its migration, by all the disc (bottom, solid line), the Hill 
sphere (dashed line b — 1 .0), the inner 80% of the Hill sphere 
(b = 0.8), 60% (b = 0.6), and 40% (b = 0.4). Bottom panel : the 
same torques as a percentage of the total torque. 

0.5 rn of the planet contains most of the mass and yet it exerts 
a negligible torque on the planet, whereas the region between 
0.5 rn and r# has a negligible mass but it exerts a large torque 
on the planet. Such a large torque cannot be permanently ex- 
erted by material bound to the planet, otherwise this material 
would not migrate at the same rate as the planet's, and it would 
ultimately escape from the planet's potential well. This suggests 
that the CPD is located inside 0.5 r#, which can be checked with 
a streamline analysis, as shown below. Alternatively, one can de- 
termine the shape and the size of the CPD with the following 
energetic approach. We denote by E tot the total specific energy 
of the fluid elements perturbed by the planet. It is the sum of the 
potential, kinetic, and thermal perturbed energies by unit mass 
of the fluid elements : 



E tot = -GM p ls' + |v - v//2 + (T- r„)/(r - 1) , 



(11) 



where v - v p is the velocity with respect to the planet, and T - 
Tq is the perturbed temperature with respect to the initial state. 
The minimum energy that the fluid elements must have to leave 
the planet's potential well with zero velocity with respect to the 
planet defines the escape energy, denoted by E esc . It reads 



£esc = -GM p ld' +(T- 7oV/(7 - 1) , 



(12) 



with d' 



4 



+ e 2 the smoothed separation between the 

planet and one of the two hyperbolic stagnation points (where 
the velocity of the flow with respect to the planet cancels out). 
Fluid elements such that E tot < E esc should be bound to the 
planet, and thus be part of the CPD, whereas those having 
E tot > E esc should orbit the central star. Said differently, the 
shape and the size of the CPD should correspond to the loca- 
tions in the disc where E tot = E esc . 

We first compare both approaches to determine the CPD 
shape with the middle-resolution run for the Saturn mass planet. 
In the top panel of Fig. [18] we depict the gas surface density at 
300 orbits, just before the planet release. The latter is located 
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Fig. 18. Top : surface density contour at 300 orbits for the 
middle-resolution run with the Saturn-mass planet. The planet 
and its Hill sphere are respectively represented by the + symbol 
and the dashed circle. The grid mesh is superimposed, as well as 
streamlines in the planet's frame. Bottom : same as the top panel, 
except the quantity (E tot - £esc)/|£escl is now displayed. The red 
solid curve shows the contour E Xoi = E esc , and the x symbol 
shows the location of the stagnation point used to evaluate the 
escape energy £e SC . 



at r = r p and <p = tp p , and its position is highlighted with a + 
symbol. The background vertical and horizontal lines show the 
grid mesh. The dashed circle represents the planet's Hill sphere. 
Streamlines (solid curves) are also superimposed to appreciate 
the shape and the size of the CPD. The streamline analysis con- 
firms that most of the mass inside the planet's Hill sphere is con- 
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fined to the CPD, and it shows that the latter has an elliptical 
shape, with a semi-major axis approximately twice as large as 
the semi-minor axis. The semi-major axis is ~ 0.6 This justi- 
fies that the material inside ~ 0.5 r# exerts a small torque on the 
planet compared to the material located between 0.5 r# and r#. 

In the bottom panel of Fig. [18] we display the quantity 
(£ t ot - £esc)/|£esd at the same time. The same streamlines as 
in the top panel are depicted. The red solid line shows the con- 
tour E tot = £e SC , the quantity E esc being calculated with Eq. ( [T2l >. 
The location of the stagnation point, depicted with a x symbol 
in this panel (at r « 0.97r p , ip x (p p ), is inferred from the stream- 
line analysis. The contour E tot = E esc can be approximated as a 
circle of radius 0.6 r#. Its location agrees with that of the CPD 
determined with the streamline analysis. 

Figure [19] is the same as Fig. [18] but for the very high- 
resolution run just before the release of the Saturn mass planet 
at 475 orbits. A close comparison between both figures reveals 
significant differences in the flow topology near the planet. In 
particular, the CPD is now very close to a circle of radius 0.5 
and its shape and its size are in very good agreement with those 
of the contour E lot = E tic . The flow also tends to show more 
complexity at higher resolution, judging from the presence of a 
single vortex slightly inside the planet's orbit. 

The C PD is slightly larger tha n in previous locally isothermal 
studies bv lD'Angelo et alJd2005l) . This probably comes from the 
structure inside the Hill sphere changeing when one changes the 
equation of state, as the pressure support is modified. With the 
energy equation, the collapse of the CPD is limited by the heat- 
ing due to adiabatic compression, which may give a wider but 
less massive CPD than in the locally isothermal case. Also, the 
smoothing length e influences the mass and size of the CPD : the 
smaller e, the deeper the potential well of the planet, and then 
the larger the CPD. The choice of the equation of state, opacity, 
and smoothing length certainly influences the size and mass of 
the CPD. The analysis performed in this section concerns only 
the specific cases that we studied in this paper, in order to com- 
pare the mass and torque repartition in the Hill sphere with the 
migration rates observed. 

7. Conclusion 

In this study, we have performed experiments to analyse the role 
of the circum-planetary disc and the gas in the Hill sphere in 
the numerical simulations of planetary migration of giant plan- 
ets. The way the material inside the Hill sphere is taken into 
account has a strong influence on the outcome of simulations in 
which the gas selfgravity is not computed. Whether or not a part 
of the Hill sphere is excluded in the computation of the force 
felt by the planet, the migration rate can vary by a factor 2 in 
type II migration, and even 3.33 in case of type III migration. 
This numerical issue is thus critical. Therefore, we would like to 
encourage authors to specify how they deal with the Hill sphere 
of giant planets in their simulations, so that their experiments 
can be reproducible by others. 

In Sect. 12.41 we analysed the problem and proposed several 
solutions. To deliver the migrating planet from the ball and chain 
effect of the CPD, it is possible to exclude a part of the Hill 
sphere from the computation of the force exerted by the disc on 
the planet. It is also possible to artificially give the CPD the ac- 
celeration felt by the planet. Both methods are not equivalent. 
In order to take the CPD into account in the perturbation of the 
PPD, the mass of the CPD can be added to the perturbing mass 
of the planet. In the second case studied here, this leads to an ac- 
celeration of the migration proportional to the mass of the CPD. 
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Fig. 19. Same as Fig. [18] but for the very-high-resolution run, 
just before the release of the planet at 475 orbits. To improve 
legibility, the grid mesh is not displayed. 

To help the planet to pull the CPD, the gravitational mass of the 
planet can be increased by Mcpd, which also accelerates the mi- 
gration as expected. These four tricks can be combined. 

We would like to point out again that the CPD is not the en- 
tire Hill sphere. As demonstrated in Sect. 16.31 its size is only 
about a half of the Hill radius only. As a consequence, exclud- 
ing all the Hill sphere of the planet is inappropriate ; in fact, our 
simulations show that this method gives the most divergent re- 
sults with respect to our reference case (full gas selfgravity com- 
puted). Unfortunately, no method matches the full selfgravity 
case. The role of the selfgravity is more complex than changing 
the inertial, perturbing, and gravitational masses of the planet. 
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The present study was conducted for a specific equation of 
state and for a specific value of the softening length. Under these 
circumstances, it seems appropriate to remove about half of the 
Hill sphere in the computation of the force felt by the planet 
from the disc, but not more than 0.6r//. Authors who consider 
different prescriptions should undertake a prior analysis of the 
CPD size to infer how they should scale this fraction in their 
case. Alternatively, accelerating the material within ~ 0.5 r# of 
the planet also provides satisfactory results, but the acceleration 
imposed on the CPD should be the one provided to the planet by 
the disc outside 0.5r H . The addition of Mqpd to the perturbing 
mass of the planet also seems to be a good idea and leads to an 
acceleration of the migration. Of course, the most reliable solu- 
tion is to compute the full gas selfgravity, with the highest res- 
olution possible. In any case, one should be aware of the strong 
influence of the CPD on planetary migration. 

Acknowledgements. A. Crida acknowledges the support through the German 
Research Foundation (DFG) grant KL 650/7. The high-resolution computations 
were performed on the hpc-bw and hpc-uni clusters of the Rechenzentrum of 
the University of Tubingen. 



References 

Baruteau, C. & Masset, F. 2008a, ApJ, 672, 1054 

Baruteau, C. & Masset, F. 2008b, in AAS/Division of Dynamical Astronomy 
Meeting, Vol. 39, AAS/Division of Dynamical Astronomy Meeting, #03.03- 
+ 

— . 2008c, ApJ, 678, 483 

Crida, A. 2006, PhD thesis, Observatoire de la Cote d'Azur, NICE, France 

Crida, A. & Morbidelli, A. 2007, MNRAS, 377, 1324 

Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 

Crida, A., Sandor, Z., & Kley, W. 2008, A&A, 483, 325 

D'Angelo, G., Bate, M. R., & Lubow, S. H. 2005, MNRAS, 358, 316 

De Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 

Goldreich, P. & Tremaine, S. 1979, ApJ, 233, 857 

Kley, W. & Crida, A. 2008, A&A, 487, L9 

Lin, D. N. C. & Papaloizou, J. 1979, MNRAS, 186, 799 

— . 1986a, ApJ, 307, 395 

— . 1986b, ApJ, 309, 846 

Masset, F. 2000a, A&AS, 141, 165 

Masset, F. 2000b, in ASP Conf. Ser. 219: Disks, Planetesimals, and Planets, ed. 

G. Garzon, C. Eiroa, D. de Winter, & T. J. Mahoney, 75-80 
Masset, F. 2006, in Tidal interactions in composite systems, ed. M.-J. Goupil & 

J.-P. Zahn 
Masset, F. S. 2001, ApJ, 558, 453 

Masset, F. S., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730 

Masset, F. S. & Papaloizou, J. C. B. 2003, ApJ, 588, 494 

Meyer- Vernet, N. & Sicardy, B. 1987, Icarus, 69, 157 

Nelson, A. F. & Benz, W. 2003, ApJ, 589, 556 

Paardekooper, S.-J. & Mellema, G. 2006, A&A, 459, L17 

Peplihski, A., Artymowicz, P., & Mellema, G. 2008, MNRAS, 386, 164 

Pierens, A. & Hure, J.-M. 2005, A&A, 433, L37 

Pierens, A. & Nelson, R. P. 2008, A&A, 482, 333 

Ward, W. R. 1986, Icarus, 67, 164 

Ward, W. R. 1991, in Lunar and Planetary Institute Conference Abstracts, 1463- 
1464 

— . 1997, Icarus, 126, 261 



